home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 011 / guard.cq / guard.c
C/C++ Source or Header  |  1985-06-16  |  2KB  |  89 lines
  1. /* guard.c
  2. *
  3. *  from the GUARD programs in Richard Karpinski's article on PARANOIA,
  4. *  DDJ, February 1985;
  5. *
  6. *  C translation by R. E. Sawyer
  7. */
  8.  
  9. #include <math.h>
  10.  
  11. main()
  12.     {
  13.     double fabs();
  14.     double               /*floating-point constants*/
  15.         one = 1.,
  16.         half = .5,
  17.         zero = 0.,
  18.         minusOne = -1.;
  19.     double               
  20.         radix,           /*calculated floating-point radix*/
  21.         precision,       /*significant digits in base radix*/
  22.         width,           /*radix^precision*/
  23.         wide,            /*first estimate of width*/
  24.         ulpOne,          /*unit in last place of just less than one*/
  25.         ulpRadix,        /*radix * ulpOne*/
  26.         oneMinus,        /*one - ulpOne  (calculated with care)*/
  27.         radixMinus,      /*radix - ulpRadix*/
  28.         s, t, u, x, y, z;/*working variables*/
  29.  
  30. /*find wide so big that adding one does not change it by one*/
  31.  
  32.     wide = one;
  33.     do
  34.         {
  35.         wide += wide;     
  36.         x = wide + one;
  37.         y = x - wide;
  38.         z = y - one;
  39.         } while (minusOne + fabs(z) < zero);
  40.  
  41. /*find the radix (or number base) as the minimum increase in wide
  42. * -- since wide is just large enough that the units place is not
  43. *    represented, a one in the last represented place is exactly
  44. *    the radix
  45. */
  46.     y = one;
  47.     do
  48.         {
  49.         radix = wide + y;
  50.         y += y;
  51.         radix -= wide;
  52.         } while (radix == zero);
  53.     printf("\n radix = %.15e", radix);
  54.  
  55. /*find precision in radix digits*/
  56.  
  57.     precision = zero;
  58.     width = one;
  59.     do
  60.         {
  61.         ++precision;
  62.         width *= radix;
  63.         y = width + one;
  64.         } while (y - width == one);
  65.     printf("\n precision = %.15e", precision);
  66.     printf("\n width = %.15e", width);
  67.     ulpOne = one / width;
  68.     printf("\n closest relative separation found is ulpOne = %.15e",
  69.             ulpOne);
  70.     oneMinus = (half - ulpOne) + half;
  71.     ulpRadix = radix * ulpOne;
  72.     radixMinus = radix - one;
  73.     radixMinus = (radixMinus - ulpRadix) + one;
  74.     x = one - ulpOne;
  75.     y = one - oneMinus;
  76.     z = one - x;
  77.     s = radix - ulpRadix;
  78.     t = radix - radixMinus;
  79.     u = radix - s;
  80.     if ((y == ulpOne) && (z == ulpOne)
  81.         && (t == ulpRadix) && (u == ulpRadix))
  82.         printf("\n add/subtract has an appropriate guard digit");
  83.     else
  84.         printf("\n add/subtract lacks a guard digit");
  85.     printf("\n");
  86.     }
  87.     
  88.     
  89.